---
title: "Replicate Table D4: Polynomial Order = 6"
author: "Arthur Zeyang Yu"
date: '2023-04-20'
output: html_document
---

```{r setup, include = TRUE}
knitr::opts_chunk$set(echo = TRUE)

# clear memory
rm(list = ls())

# set working directory in Thinkpad
setwd(getwd())

# load packages
library(MASS)
library(Matrix)
library(data.table)
library(stats)
library(plm)
library(xtable)
library(slam)
library(gurobi)

# run non linear programming R code
source("r_ate_at_code_nonlp.R")
```

# prepare: Bernstein polynomials function
```{r}
# compute: the integral of Bernstein polynomials
integral <- function(k,K){
  
  # define bernstein polynomials
  bernstein <- function(u){
    b <- choose(K,k) * u^k * (1-u)^(K-k)
    return(b)
  }
  
  i1 <- matrix(0,1,2)
  
  i1[1,1] <- integrate(bernstein, 0, p0)[[1]]
  i1[1,2] <- integrate(bernstein, 0, p1)[[1]]
  
  i0 <- matrix(0,1,2)
  
  i0[1,1] <- integrate(bernstein, p0, 1)[[1]]
  i0[1,2] <- integrate(bernstein, p1, 1)[[1]]
  
  i <- cbind(i1,i0)
  
  return(i)
  
}

# integral of Bernstein for each k,d,z when basis is 4th degreee
K0 <- 6
k0 <- t(seq(0,K0))
i0 <- mapply(integral, k0, K0)

id1 <- i0[1:2,]
id0 <- i0[3:4,]
```

# construct constraints in the linear program
```{r}
# construct gamma_iv: constraint
## rmk1: the full weights for gamma_iv is in mst ecta paper eq (9)
## rmk2: taking expectation over Z, with value = 0, 1
gamma_iv1 <- (0.5 * (c(0, 1) - e_z)/cov_dz) %*% id1
gamma_iv0 <- (0.5 * (c(0, 1) - e_z)/cov_dz) %*% id0
gamma_iv <- cbind(gamma_iv0, gamma_iv1)

# construct gamma_ols: constraint
## rmk1: the full weights for gamma_iv is in mst ecta paper eq (9)
## rmk2: taking expectation over Z, with value = 0, 1
gamma_ols1 <- (0.5 * (c(0, 1) - e_d)/var_d) %*% id1
gamma_ols0 <- (0.5 * (c(0, 1) - e_d)/var_d) %*% id0
gamma_ols <- cbind(gamma_ols0, gamma_ols1)

# construct gamma_e_ydz: constraint
## rmk1: there are four constraints given we have binary iv and treatment
gamma_eyd1z1_1 <- matrix(0, 1, 7)
gamma_eyd1z1_0 <- (c(0, 1) %*% id1) * e_z
gamma_eyd1z1 <- cbind(gamma_eyd1z1_1, gamma_eyd1z1_0)

gamma_eyd1z0_1 <- matrix(0, 1, 7)
gamma_eyd1z0_0 <- (c(1, 0) %*% id1) * (1 - e_z)
gamma_eyd1z0 <- cbind(gamma_eyd1z0_1, gamma_eyd1z0_0)

gamma_eyd0z1_1 <- (c(0, 1) %*% id0) * e_z
gamma_eyd0z1_0 <- matrix(0, 1, 7)
gamma_eyd0z1 <- cbind(gamma_eyd0z1_1, gamma_eyd0z1_0)

gamma_eyd0z0_1 <- (c(1, 0) %*% id0) * (1 - e_z)
gamma_eyd0z0_0 <- matrix(0, 1, 7)
gamma_eyd0z0 <- cbind(gamma_eyd0z0_1, gamma_eyd0z0_0)

# define: beta_tsls, note that in present set up, wald = iv = tsls
beta_tsls <- beta_iv

# compute: beta_ols
e_yd <- e_y1_d1 * e_d
beta_ols <- (e_yd - e_y * e_d)/var_d
```

# construct weights in the objective function
```{r}
# compute: weights for E[Y(0) | D(1) = D(0) = 1] (always-taker)
gamma_ate_at0 <- matrix(0, 1, 7)
gamma_ate_at1 <- (c(1, 0) %*% id1) / p0
gamma_ate_at <- cbind(gamma_ate_at1, gamma_ate_at0)

# compute: weights for E[Y(1) | D(1) = D(0) = 0] (never-taker)
gamma_ate_nt0 <- (c(0, 1) %*% id0) / (1 - p1)
gamma_ate_nt1 <- matrix(0, 1, 7)
gamma_ate_nt <- cbind(gamma_ate_nt1, gamma_ate_nt0)
```

# produce partial identification results
```{r}
# Row 5: bounds with iv as constraint
model1 <- list()
model1$A <- rbind(gamma_iv,Diagonal(14),Diagonal(14))
model1$obj <- gamma_ate_at
model1$modelsense <- 'max'
model1$rhs <- c(beta_iv,matrix(0,14,1),matrix(1,14,1))
model1$sense <- c('=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

model1 <- list()
model1$A <- rbind(gamma_iv,Diagonal(14),Diagonal(14))
model1$obj <- gamma_ate_at
model1$modelsense <- 'min'
model1$rhs <- c(beta_iv,matrix(0,14,1),matrix(1,14,1))
model1$sense <- c('=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

model1 <- list()
model1$A <- rbind(gamma_iv,Diagonal(14),Diagonal(14))
model1$obj <- gamma_ate_nt
model1$modelsense <- 'max'
model1$rhs <- c(beta_iv,matrix(0,14,1),matrix(1,14,1))
model1$sense <- c('=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

model1 <- list()
model1$A <- rbind(gamma_iv,Diagonal(14),Diagonal(14))
model1$obj <- gamma_ate_nt
model1$modelsense <- 'min'
model1$rhs <- c(beta_iv,matrix(0,14,1),matrix(1,14,1))
model1$sense <- c('=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

# Row 8: bounds with iv + ols in constraint
model1 <- list()
model1$A <- rbind(gamma_iv, gamma_ols, Diagonal(14),Diagonal(14))
model1$obj <- gamma_ate_at
model1$modelsense <- 'max'
model1$rhs <- c(beta_iv, beta_ols, matrix(0,14,1),matrix(1,14,1))
model1$sense <- c('=','=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

model1 <- list()
model1$A <- rbind(gamma_iv, gamma_ols, Diagonal(14),Diagonal(14))
model1$obj <- gamma_ate_at
model1$modelsense <- 'min'
model1$rhs <- c(beta_iv, beta_ols, matrix(0,14,1),matrix(1,14,1))
model1$sense <- c('=','=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

model1 <- list()
model1$A <- rbind(gamma_iv, gamma_ols, Diagonal(14),Diagonal(14))
model1$obj <- gamma_ate_nt
model1$modelsense <- 'max'
model1$rhs <- c(beta_iv, beta_ols, matrix(0,14,1),matrix(1,14,1))
model1$sense <- c('=','=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

model1 <- list()
model1$A <- rbind(gamma_iv, gamma_ols, Diagonal(14),Diagonal(14))
model1$obj <- gamma_ate_nt
model1$modelsense <- 'min'
model1$rhs <- c(beta_iv, beta_ols, matrix(0,14,1),matrix(1,14,1))
model1$sense <- c('=','=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

# Row 11: bounds for AT by add E[1{D = 0, Z = 0} Y] in constraint
model1 <- list()
model1$A <- rbind(gamma_eyd0z0, Diagonal(14), Diagonal(14))
model1$obj <- gamma_ate_at
model1$modelsense <- 'max'
model1$rhs <- c(e_d0z0, matrix(0,14,1), matrix(1,14,1))
model1$sense <- c('=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

model1 <- list()
model1$A <- rbind(gamma_eyd0z0, Diagonal(14), Diagonal(14))
model1$obj <- gamma_ate_at
model1$modelsense <- 'min'
model1$rhs <- c(e_d0z0, matrix(0,14,1), matrix(1,14,1))
model1$sense <- c('=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

model1 <- list()
model1$A <- rbind(gamma_eyd0z0, Diagonal(14), Diagonal(14))
model1$obj <- gamma_ate_nt
model1$modelsense <- 'max'
model1$rhs <- c(e_d0z0, matrix(0,14,1), matrix(1,14,1))
model1$sense <- c('=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

model1 <- list()
model1$A <- rbind(gamma_eyd0z0, Diagonal(14), Diagonal(14))
model1$obj <- gamma_ate_nt
model1$modelsense <- 'min'
model1$rhs <- c(e_d0z0, matrix(0,14,1), matrix(1,14,1))
model1$sense <- c('=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

# Row 14: bounds for AT by add E[1{D = 1, Z = 0} Y] in constraint
model1 <- list()
model1$A <- rbind(gamma_eyd1z0, Diagonal(14), Diagonal(14))
model1$obj <- gamma_ate_at
model1$modelsense <- 'max'
model1$rhs <- c(e_d1z0, matrix(0,14,1), matrix(1,14,1))
model1$sense <- c('=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

model1 <- list()
model1$A <- rbind(gamma_eyd1z0, Diagonal(14), Diagonal(14))
model1$obj <- gamma_ate_at
model1$modelsense <- 'min'
model1$rhs <- c(e_d1z0, matrix(0,14,1), matrix(1,14,1))
model1$sense <- c('=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

model1 <- list()
model1$A <- rbind(gamma_eyd1z0, Diagonal(14), Diagonal(14))
model1$obj <- gamma_ate_nt
model1$modelsense <- 'max'
model1$rhs <- c(e_d1z0, matrix(0,14,1), matrix(1,14,1))
model1$sense <- c('=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

model1 <- list()
model1$A <- rbind(gamma_eyd1z0, Diagonal(14), Diagonal(14))
model1$obj <- gamma_ate_nt
model1$modelsense <- 'min'
model1$rhs <- c(e_d1z0, matrix(0,14,1), matrix(1,14,1))
model1$sense <- c('=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

# Row 17: 2 moments (E[1{D = 0, Z = z} Y] w\ z \in {0, 1}) in constraint moments
model1 <- list()
model1$A <- rbind(gamma_eyd0z0, gamma_eyd0z1, Diagonal(14),Diagonal(14))
model1$obj <- gamma_ate_at
model1$modelsense <- 'max'
model1$rhs <- c(e_d0z0, e_d0z1, matrix(0,14,1),matrix(1,14,1))
model1$sense <- c('=','=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

model1 <- list()
model1$A <- rbind(gamma_eyd0z0, gamma_eyd0z1, Diagonal(14),Diagonal(14))
model1$obj <- gamma_ate_at
model1$modelsense <- 'min'
model1$rhs <- c(e_d0z0, e_d0z1, matrix(0,14,1),matrix(1,14,1))
model1$sense <- c('=','=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

model1 <- list()
model1$A <- rbind(gamma_eyd0z0, gamma_eyd0z1, Diagonal(14),Diagonal(14))
model1$obj <- gamma_ate_nt
model1$modelsense <- 'max'
model1$rhs <- c(e_d0z0, e_d0z1, matrix(0,14,1),matrix(1,14,1))
model1$sense <- c('=','=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

model1 <- list()
model1$A <- rbind(gamma_eyd0z0, gamma_eyd0z1, Diagonal(14),Diagonal(14))
model1$obj <- gamma_ate_nt
model1$modelsense <- 'min'
model1$rhs <- c(e_d0z0, e_d0z1, matrix(0,14,1),matrix(1,14,1))
model1$sense <- c('=','=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

# Row 20: 2 moments (E[1{D = 1, Z = z} Y] w\ z \in {0, 1}) in constraint moments
model1 <- list()
model1$A <- rbind(gamma_eyd1z0, gamma_eyd1z1, Diagonal(14),Diagonal(14))
model1$obj <- gamma_ate_at
model1$modelsense <- 'max'
model1$rhs <- c(e_d1z0, e_d1z1, matrix(0,14,1),matrix(1,14,1))
model1$sense <- c('=','=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

model1 <- list()
model1$A <- rbind(gamma_eyd1z0, gamma_eyd1z1, Diagonal(14),Diagonal(14))
model1$obj <- gamma_ate_at
model1$modelsense <- 'min'
model1$rhs <- c(e_d1z0, e_d1z1, matrix(0,14,1),matrix(1,14,1))
model1$sense <- c('=','=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

model1 <- list()
model1$A <- rbind(gamma_eyd1z0, gamma_eyd1z1, Diagonal(14),Diagonal(14))
model1$obj <- gamma_ate_nt
model1$modelsense <- 'max'
model1$rhs <- c(e_d1z0, e_d1z1, matrix(0,14,1),matrix(1,14,1))
model1$sense <- c('=','=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval

model1 <- list()
model1$A <- rbind(gamma_eyd1z0, gamma_eyd1z1, Diagonal(14),Diagonal(14))
model1$obj <- gamma_ate_nt
model1$modelsense <- 'min'
model1$rhs <- c(e_d1z0, e_d1z1, matrix(0,14,1),matrix(1,14,1))
model1$sense <- c('=','=',
                  '>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=','>=',
                  '<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=','<=')
result1 <- gurobi(model1)
result1$objval
```
